clear all; clc; close all;

tstart = tic % время вычисления программы

L = 100e-6;
C = 33e-12;
R = 120;

w0 = 1/sqrt(L*C);
dw = w0/100;
wmin = w0 - w0*0.5;
wmax = w0 + w0*0.5;
w = wmin:dw:wmax;

a = [L R 1/C];% [a2 a1 a0]
b = [1 0];% [b1 b0]
% Находим ПФ, строим АЧХ и ФЧХ
figure(1)
freqs(b,a,w); % Вывести АЧХ/ФЧХ, приняв по оси частот w
H = freqs(b,a,w); % Посчитать ПФ, но график не выводить

figure(2)
subplot(2,1,1)
semilogy(w/2/pi/1e6, abs(H));
xlabel('f, MHz')
ylabel('|H|')
% 
subplot(2,1,2)
plot(w/2/pi/1e6, unwrap(angle(H)));
xlabel('f, MHz')
ylabel('arg H')
% Находим ИХ
sys = tf(b,a);
[y,t] = impulse(sys);
figure(3)
plot(t, y);
xlabel('t, sec')
ylabel('h(t)')

% Находим цифровой аналог
Fs = 10*wmax/2/pi;
[bz, az] = bilinear(b,a,Fs); % Получаем цифровые коэффициенты
[Hz, wz] = freqz(bz, az);
figure;
plot(w/2/pi/1e6, abs(H), Fs*wz/2/pi/1e6, abs(Hz));
xlabel('f, MHz')
ylabel('|H|')
grid on

%формируем шум
 
Td = 1/ Fs;
Tmod = 1000*Td;
t = 0:Td:Tmod; 

stdn = 13; 
n = stdn*randn(1, length(t));
y = filter(bz,az,n);

figure
plot(t,[n;y]) % упрощение plot(t, n, t, y)
xlabel('t, sec')
ylabel('n(t),y(t)')

%докажем,что нормальный
%hist(n,100)

%спектр шума на вх и вых


nf = fft(n);
yf = fft(y);

mn = sqrt(mean (abs(nf).^2))*2.5;
my = sqrt(mean (abs(yf).^2))*10;
nf = nf/mn;
yf = yf/my;

f = 0:1/Tmod:1/Td;
figure;
plot(f/1e6, [abs(nf); abs(yf)], w/2/pi/1e6, abs(H)/max(abs(H)), 'LineWidth',2);
xlim([0 Fs/1e6/2])
grid on
xlabel('f, MHz')

% время вычисления программы
dt = toc(tstart);
fprintf('dt = %f ms\n', dt*1e3);



